This notebook is Part 3 of an exploration to visualize results of hydrologic models. In Part 1 we built custom HTML widgets using D3.js, and in Part 2 we looked at rendering water fluxes using Sankey diagrams. Here we test multiple libraries to generate hillshade (3D) views of river basins and water infrastructure, in particular we want to compare leaflet, Three.js and WebGL implementations.
Aside from rendering topography and water streams in 3D (and potentially other covariate layers), our objective is to overlay custom labels to illustrate the water cycle.
Another objective is to provide basin and sub-basin statistics, starting with precipitation, ET, soil moisture, and adding other covariates, such as population, land use allocation, and crop allocation.
Some inspiration below:
- Tylor Morgan Tutorial: Adding OSM Data to Rayshader Maps
- Flood depth on 3D terrain
- NASA Web WorldWind
Show code
list.files("./fig", full.names=TRUE)[c(2,7)] %>%
knitr::include_graphics()
WebGL: Sample Scene Processing with Rayshader1
Data Acquisition
We’ll experiment with a smaller sample scene of the Selingue Dam in the Niger River basin.
basin <- shapefile("~/Projects/2021-iwmi/data/mli/srtm/mli_basin.shp")
zoi <- shapefile("~/Projects/2021-iwmi/data/mli/srtm/zoi.shp")
ext <- extent(zoi)
center <- coordinates(zoi)
# Get CGIAR SRTM DEM at 90m
srtm <- getData("SRTM", lon=center[,1], lat=center[,2], path="_data") %>%
crop(zoi) %>%
mask(zoi)
srtm
class : RasterLayer
dimensions : 838, 973, 815374 (nrow, ncol, ncell)
resolution : 0.0008333333, 0.0008333333 (x, y)
extent : -8.6725, -7.861667, 11.06917, 11.7675 (xmin, xmax, ymin, ymax)
crs : +proj=longlat +datum=WGS84 +no_defs
source : memory
names : srtm_35_10
values : 326, 476 (min, max)The basin covers a large area, so we would need 8 SRTM tiles, but 1 is enough for a proof of concept. Next we’ll get a satellite basemap.
plot(terra::vect(basin), col=pal[2], border=pal[1], lwd=2,
main="Niger River Basin (Mali)")
plot(zoi, lty=3, col=alpha(pal["red"], .6), border=pal["red"], lwd=2, add=T)
text(-8, 10.5, "Selingue Dam\n(Mali)", col=pal["red"], cex=.7, font=2)
grid()
plot(terra::ext(zoi),
main="Selingue Dam (Mali) - ESRI World Imagery")
plotRGB(bmap, add=T)
grid(col="white")
plot(terra::rast(srtm),
main="Selingue Dam (Mali) - SRTM 90m")
grid(col="white")
Scene Rendering
Next we convert the 2 rasters to a matrix format that’s compatible with Rayshader hillshading methods.
Show code
# Convert rasters to rayshader matrix format
srtm_array <- raster_to_matrix(srtm)
# Convert sat basemap to matrix (test)
r <- raster_to_matrix(bmap$red)
g <- raster_to_matrix(bmap$green)
b <- raster_to_matrix(bmap$blue)
bmap_array <- array(0, dim=c(nrow(r), ncol(r), 3))
bmap_array[,,1] <- r/255
bmap_array[,,2] <- g/255
bmap_array[,,3] <- b/255
bmap_array %>%
aperm(c(2,1,3)) %>%
# Stretch contrast
rescale(to=c(0,1)) -> bmap_array
srtm_water <- srtm_array
srtm_water[srtm_water < 353] <- 0
basemap_sat <- srtm_array %>%
height_shade() %>%
add_overlay(bmap_array) %>%
add_shadow(ray_shade(srtm_array, zscale=90)) %>%
add_water(detect_water(srtm_water), color=alpha(pal["blue"], 0.4))
basemap_sat %>% plot_map()
Figure 1: Hillshade Basemap of Selingue Dam (Niger River basin)
That doesn’t look very clear, so instead we’ll create a basemap, not using the satellite image but a built-in desert texture.
base_map <- srtm_array %>%
height_shade() %>%
add_overlay(sphere_shade(srtm_array, texture="desert",
zscale=4, colorintensity=5), alphalayer=0.5) %>%
add_shadow(lamb_shade(srtm_array, zscale=6), 0) %>%
add_shadow(ambient_shade(srtm_array), 0) %>%
add_shadow(texture_shade(srtm_array, detail=8/10, contrast=9, brightness=11), 0.1) %>%
add_water(detect_water(srtm_water), color=alpha(pal["blue"], 0.4))
saveRDS(base_map, "./_data/base_map.rds")
Show code
Figure 2: Step 1 - Hillshade Basemap of Selingue Dam
Looks better, so let’s acquire and overlay spatial features from OSM. Note that water-related features for Mali are scarcely available from OSM, but network streams are recorded.
osm_roads <- opq(bbox(zoi)) %>%
add_osm_feature("highway") %>%
osmdata_sf()
osm_water = opq(bbox(zoi)) %>%
add_osm_feature("waterway") %>%
osmdata_sf()
osm_place = opq(bbox(zoi)) %>%
add_osm_feature("place") %>%
osmdata_sf()
road_layer <- generate_line_overlay(
dplyr::filter(osm_roads$osm_lines, highway %in% c("primary", "secondary")),
extent=ext, srtm_array, linewidth=5, color=pal["black"])
water_layer <- generate_line_overlay(
osm_water$osm_lines,
extent=ext, srtm_array, linewidth=3, color=pal["blue"])
place_layer <- generate_label_overlay(
dplyr::filter(osm_place$osm_points, !is.na(name) & nchar(name)<10),
extent=ext, heightmap=srtm_array,
font=2, text_size=1.6, point_size=1.6, color=pal["black"],
halo_color="white", halo_expand=2, halo_blur=1, halo_alpha=.9, seed=1,
data_label_column="name")
scene <- base_map %>%
#scene <- basemap_sat %>%
add_overlay(road_layer) %>%
add_overlay(water_layer, alphalayer=1) %>%
add_overlay(place_layer)
saveRDS(scene, "./_data/scene.rds")
Finally we’ll use WebGL to render this scene in 3D. We also test adding polygon (bar) annotations.
amb_layer <- ambient_shade(srtm_array, zscale=1/5)
scene2 <- srtm_array %>%
height_shade() %>%
add_shadow(texture_shade(srtm_array, detail=8/10, contrast=9, brightness=11), 0) %>%
add_shadow(amb_layer, 0) %>%
add_overlay(road_layer) %>%
add_overlay(water_layer, alphalayer=1) %>%
add_overlay(place_layer)
saveRDS(scene2, "./_data/scene2.rds")
Show code
scene <- readRDS("./_data/scene.rds")
scene2 <- readRDS("./_data/scene2.rds")
scene %>% plot_3d(alt, zscale=20,
theta=30, phi=20, fov=0, zoom=0.5)
# Add polygon annotations
xyz <- sf::read_sf("~/Projects/2021-iwmi/data/mli/srtm/xyz.shp")
render_polygons(xyz, ext, data_column_top="z",
scale_data=1, color=alpha(pal["orange"], 0.8),
lit=F, light_intensity=0.01, clear_previous=T)
rglwidget()
Figure 3: Step 2: Interactive 3D Scene of Selingue Dam
The generated scene (PNG image and associated JS code) is 3.5MB.
WebGL: Full Scene with Rayshader
Following the approach above, we experiment with a 3D scene for the entire Niger River basin area used in the analysis.
Data Acquisition
# Basin boundaries
zoi <- shapefile("~/Projects/2021-iwmi/data/mli/srtm/mli_basin.shp")
ext <- extent(zoi)
center <- coordinates(zoi)
# Admin2 boundaries
adm <- sf::st_read(
"~/Projects/2021-iwmi/data/mli/srtm/mli_adm2_lines.shp")
places <- sf::st_read(
"~/Projects/2021-iwmi/data/mli/srtm/mli_adm2_centroids.shp")
# Get CGIAR SRTM DEM (upsampled)
alt <- lapply(c("MLI", "GIN", "CIV"),
function(x) getData("alt", country=x, mask=FALSE) %>%
crop(zoi) %>%
mask(zoi)
)
alt <- mosaic(alt[[1]], alt[[2]], alt[[3]], fun=mean)
alt <- mask(crop(alt, zoi), zoi)
writeRaster(alt, "~/Projects/2021-iwmi/data/mli/srtm/mli_alt.tif",
overwrite=T)
alt
# Get ESA WorldCover (clipped in QGIS)
luc <- raster("~/Projects/2021-iwmi/data/mli/srtm/ESA WorldCover clip.tiff")
# Note that I can't find a colormap for this raster, so switching to ESA CCI instead
luc <- raster("~/Projects/2021-iwmi/data/mli/srtm/mli_esa_cci_300m.tif")
# And its colormap
pal.luc <- fread("~/Maps/ESA/ESA CCI Colormap.txt")
setnames(pal.luc, c("value", "R", "G", "B", "A", "label"))
pal.luc[, hex := rgb(R, G, B, maxColorValue=255)]
# OSM basemap
bmap <- maptiles::get_tiles(terra::ext(zoi), "OpenStreetMap.HOT", zoom=6) %>%
stack() %>%
crop(zoi) %>%
mask(zoi)
# We'll use GloRIC v10 stream network instead of OSM
# Features are filtered to 'Class_hydr > 12'
gloric <- sf::st_read("~/Projects/2021-iwmi/data/mli/srtm/mli_gloric_filtered.shp")
save(zoi, ext, adm, places, alt, luc, pal.luc, bmap, gloric,
file="./_data/osm.RData")
load("./_data/osm.RData")
plot(terra::rast(alt),
main="Niger River Basin (Mali) - SRTM 1km")
plot(adm, lty=3, col=pal["black"], lwd=1, add=T)
plot(gloric, col=pal["blue"], lwd=.4, add=T)
grid()
plot(terra::rast(luc), col=pal.luc$hex,
breaks=pal.luc$value, legend=F,
main="Niger River Basin (Mali) - ESA CCI 300m")
plot(adm, lty=3, col=pal["black"], lwd=1, add=T)
grid()
Scene Rendering
Same configuration as above.
Show code
# Convert rasters to rayshader matrix format
srtm_array <- raster_to_matrix(alt)
# Convert sat basemap to matrix (test)
r <- raster_to_matrix(bmap$red)
g <- raster_to_matrix(bmap$green)
b <- raster_to_matrix(bmap$blue)
bmap_array <- array(0, dim=c(nrow(r), ncol(r), 3))
bmap_array[,,1] <- r/255
bmap_array[,,2] <- g/255
bmap_array[,,3] <- b/255
bmap_array %>%
aperm(c(2,1,3)) %>%
# Stretch contrast
rescale(to=c(0,1)) -> bmap_array
base_map <- srtm_array %>%
height_shade() %>%
add_overlay(sphere_shade(srtm_array, texture="desert"), alphalayer=0.5) %>%
add_shadow(lamb_shade(srtm_array, zscale=6), 0) %>%
add_shadow(ambient_shade(srtm_array), 0)
saveRDS(base_map, "./_data/base_map_basin.rds")
Figure 4: Step 1 - Hillshade Basemap Niger River Basin
admin_layer <- generate_line_overlay(adm,
extent=ext, srtm_array, linewidth=3, lty=1, color=pal["light"])
water_layer <- generate_line_overlay(gloric,
extent=ext, srtm_array, linewidth=3, color=pal["light-blue"])
place_layer <- generate_label_overlay(places,
extent=ext, heightmap=srtm_array,
font=2, text_size=1.6, point_size=0, color=pal["light"],
halo_color=pal["black"], halo_expand=2, halo_blur=2, halo_alpha=.6, seed=1,
data_label_column="ADM2_NAME")
scene <- base_map %>%
#scene <- basemap_luc %>%
add_overlay(admin_layer) %>%
add_overlay(water_layer)
#add_overlay(place_layer)
saveRDS(scene, "./_data/scene_basin.rds")
Finally we’ll use WebGL to render this scene in 3D (place names need additional work).
Show code
scene <- readRDS("./_data/scene_basin.rds")
scene %>% plot_3d(resize_matrix(srtm_array, 1/5), zscale=80,
theta=30, phi=50, fov=0, zoom=0.5)
rglwidget()
Figure 5: Step 2: Interactive 3D Scene of Niger River Basin
Show code
3D Processing with Three.js
Similar implementation using Three.js libraries and viewer instead of WebGL. The scene is 12 MB in size (using a low resampling rate), it was generated using the Qgis2threejs plugin. The main difference is we’re showing ESA CCI land cover classification as an image overlay instead of a texture.
Show code
knitr::include_graphics("./fig/threejs_luc.png")
Figure 6: 3D Scene of Niger River Basin with Land Cover Overlay
This widget is rendered here.